# topic 11 h
#  Load the functions we will need in this script
source( "../assess_normality.R")
source( "../pop_sd.R")
#
#  Let us look at the binomial probabilities.
#  We did this earlier in Topic 11b but there
#  we used the function pbinom.  So, if we had
#  an problem where the probability of a single
#  success is 0.68, then we could ask for
#  the probability on 36 trials of getting
#  23 or fewer successes.
pbinom( 23, 36, 0.68)
#
#  Now, we will try this as an experiment
#  and we will do the experiment 10,000 times.
L1 <- 1:10000
for (i in 1:10000) {
  # do a 36 trial experiment
  how_many <- 0
  for (j in 1:36 ) {
    this_single <- sample(1:2,1,prob=c(0.68,0.32))
    if( this_single == 1 ){ 
       how_many <- how_many+1}
  }
  # now record the result as the number of successes
  L1[i] <-  how_many
}
# then let us look at this population of 
# 36 trial experiments
summary( L1 )
boxplot( L1, horizontal=TRUE )
hist( L1, breaks=seq(10,36,1), xaxp=c(10,36,13) )
assess_normality( L1 )
#
#  so this looks like a normal distribution
mean(L1 )
pop_sd(L1 )
#  let us try to use this to solve our 
#  earlier problem.  We know that the mean
#  of a binomial is n*p and the standard 
#  deviation of the number of successes is
#  sqrt(n*p*(1-p))
pnorm( 23, mean=36*0.68, 
       sd=sqrt(36*0.68*(1-0.68)))
# compare this to the earlier result
pbinom( 23, 36, 0.68)
# We can make an adjustment here because we
# really want to use 23.5
pnorm( 23.5, mean=36*0.68, 
       sd=sqrt(36*0.68*(1-0.68)))

#  But it seems silly to use this normal
#  approximation when we have pbinom() which
#  just gives the computed value.  The normal
#  approximation was useful years ago, before 
#  we had calculators and computers.

#
#  The steps in doing a normal approximation
#  for the binomial distribution are always the
#  same.  So we might as well capture those
#  steps in a function, pnormbino()

pnormbino <- function( num_success, p, n, lower.tail=TRUE)
{
  if(p*n < 10) {return("n*p < 10, will not compute this")}
  if(n*(1-p) < 10)
  {return("n*(1-p)<10, will not compute this")}
  nbsd <- sqrt( n*p*(1-p))
  prob <- pnorm(num_success, n*p, nbsd, lower.tail)
  return( prob )
}
pnormbino( 23.5,0.68,36 )

#  we should probably note that pnormbino is NOT
#  in our list of functions in our parent
#  directory so source will not work for it.
#  However, we do have pbinom() so pnormbino()
#  is not really needed.